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OO 1 The Kepler Mission is monitoring the brightness of ~ 150,000 stars searching 

o: 

for evidence of planetary transits. As part of the "Hunt for Exomoons with Ke- 
> '. pier" (HEK) project, we report a planetary system with two confirmed plan- 

•i-H ■ 

ets and one candidate planet discovered using the publicly available data for 



KOI-872. Planet b transits the host star with a period P b = 33.6 d and exhibits 
large transit timing variations indicative of a perturber. Dynamical model- 
ing uniquely detects an outer nontransiting planet c near the 5:3 resonance 
(P c = 57.0 d) of mass 0.37 times that of Jupiter. Transits of a third planetary 
candidate are also found: a 1.7-Earth radius super-Earth with a 6.8 d period. 
Our analysis indicates a system with nearly coplanar and circular orbits, rem- 
iniscent of the orderly arrangement within the solar system. 
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If a planet's orbit is viewed nearly edge-on, the planet may transit over the disk of its host star 
and periodically block a small fraction of the starlight. The planet's presence is then revealed 
by a small and repetitive decrease of the host star's brightness during transits. The transit light 
curve is characterized by the time of transit minimum r, the transit depth 5, the total duration 
T u , and the partial duration T 23 (1). A precise measurement of these terms allows an observer to 
infer the physical properties of the system, such as the radius ratio, p = Rp/R*, transit impact 
parameter, bp, and scaled semi-major axis, ap/i?*, where i?* is the star's physical radius and 
the subscript P denotes "planet". 

For a planet following a strictly Keplerian orbit, the spacing, timing and other properties of 
the transit light curve should be unchanging in time. Several effects, however, can produce de- 
viations from the Keplerian case so that the spacing of r is not strictly periodic and/or T 14 varies 
from transit to transit. Such changes are known as transit timing variations (TTVs) and transit 
duration variations (TDVs), respectively. TTVs are particularly sensitive to gravitational per- 
turbations from additional planets orbiting the host star (2-4) and distant large moons orbiting 
the transiting planet (5-7). 

As part of the "Hunt for Exomoons with Kepler" (HEK) project (8), we analyzed the pub- 
licly available Kepler data up to Quarter 6 (Q6; released on January 7, 2012). At the time of 
writing, the 33. 6-day -period planetary candidate Kepler-object-of-interest 872.01 (KOI-872.01) 
is the only known candidate in the system (9). The candidate was identified through HEK's 
target selection procedure as a high priority object because of the presence of visual transit 
anomalies and the dynamical capacity to host a moon. 

We detrended the raw Kepler photometry of KOI-872, covering 15 transits, using a harmonic 
filter and an exponential decay ramp correction (10). We tested several models to explain the 
photometry using the multimodal nested sampling algorithm MultiNest (11, 12), designed 
to compute the Bayesian evidence for each model. The favored model was found to be that of 
a planet undergoing TTVs (each transit has a unique r but common parameters for all other 
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terms; model Jv[ T ). This model is preferred over that of a planet on a linear ephemeris, M. p, at 
a confidence of 43.7 a (Fig.CQ). 

We computed TTVs relative to the linear ephemeris derived from M. P (P b = 33.60134 ± 
0.00021 d and r = 2455053.2826 ± 0.0014 BJD UTC , where BJD UTC is understood to be 
barycentric Julian date in coordinated universal time). The results (Fig. [2]) indicate that KOT 
872.01 exhibits large (2 hours) and complex TTVs with a dominant period of about ~5.6 transit 
cycles (~190d). These are among the largest TTVs ever detected. A model including TDVs is 
disfavored relative to the TTV-only model at a confidence of 17.5 a. 

Because a moon should induce TDVs and TTVs, the lack of TDVs does not favor a moon 
hypothesis. Further, Q4-Q6 photometry show complex stellar activity, which may be respon- 
sible for the initially identified visual anomalies. In support of these arguments, we found a 
planet-with-moon model (M. m) inadequate to explain the measured TTVs (Fig. |2fc). 

A distant stellar companion or secular/resonant perturbations from a planetary compan- 
ion also cannot explain KOI-872.01's TTVs, because these effects would produce sinusoidal 
patterns with a very long period (13). Moreover, other known TTV effects, such as parallax 
effects (14), the Applegate effect (15), and stellar proper motion (16), can be ruled out because 
they are unable to produce such a large TTV amplitude. 

We applied the TTV inversion method described in (17) to test whether the observed TTVs 
are consistent with short-period perturbations from a planetary companion and whether a unique 
set of parameters can be derived to describe the physical and orbital properties of that compan- 
ion. The short-period perturbations are small variations around the mean Keplerian orbit of 
a planet with characteristic periods comparable to the planet's orbital period. The inversion 
method is based on perturbation theory, which greatly speeds up the computation of the timing 
and duration of individual transits (18). 

We tested orbits with periods between 1 day and 10 years, including the cases of highly 
eccentric and/or retrograde planets (19). The identified solutions were fine tuned using a precise 
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iV-body integrator. We used the downhill simplex method to search for the minimum of \ 2 = 
Yuj=-i(fito,j ~ fitcjY/cji where n = 15 is the number of transits, Stoj and Stcj are the 
observed and calculated TTVs, and <7j is the uncertainty of 5t j. Local minima in x 2 were 
tested for physical plausibility, including a stability test where the orbits were tracked over 
10 9 yr with a symplectic integrator (20). 

All except two solutions can be ruled out because they are either dynamically unstable or 
have Xmin > 60. The best-fit solution (hereafter si) fits the data extremely well: Xminsi = 3-4 
for n — m = 15 — 7 = 8 degrees of freedom (DOF), where m = 7 parameters of the perturbing 
planet are the mass ratio M c /M*, period P c , eccentricity e c , inclination I c , nodal longitude Q c , 
pericenter longitude w c , and mean longitude A c . The inclination I is defined relative to the 
reference plane that is 90° tilted to the sky plane, and rotated so that £l b = 270° (all longitudes 
measured relative to the line of sight; hereafter, transit reference system). The orbital inclination 
of KOI-872.01 relative to the transit plane is h = 0.96°, as determined from the transit fit. 

The second best-fit solution (s2) has Xmins2 = 20.3. It corresponds to a planet with M c ~ 
1.8 x 10~ 3 M*, P c ~ 81.7 days (i.e. P c /P b = 2.43, just inside the 5:2 resonance), e c ~ 0.03 and 
I c ~ 10°. s2 can be ruled out because the goodness-of-fit for the TTVs is significantly worse 
at Ax 2 = 16.9; also, as I c is relatively large, s2 implies strong TDVs that are inconsistent with 
the observations (Fig. |2fc) (21). 

Therefore the transit variations of KOI-872.01 can only be fit by si. This was by no means 
expected or guaranteed because the short-periodic TTVs produced by the interacting planets 
represent only a very specific subset of astrophysical signals. This can be demonstrated by 
scrambling the TTV datapoints and applying the TTV inversion method to the scrambled data. 
We were unable to find a plausible planet solution for any of the attempted (thousand) trials 
with the scrambled data. This is a strong indication that KOI-872 is a real system of at least two 
planets (22) (Table Q). 

The scaled mass of KOI-872.02 inferred from TTVs of KOI-872.01 is M c /M* = 3.97 x 
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10~ 4 . With M* ~ 0.9 M Q , obtained from spectroscopy, this gives M c ~ 0.37 Mj, or ~ 1.3 
Saturn masses. The mass of KOI-872.01, M&, cannot be constrained from TTVs because the 
short-periodic TTVs are practically independent of M b (3, 18). The stability requirements imply 
that Mb < 6 Mj, because the orbits are dynamically unstable with a more massive transiting 
body. This mass limit along with our vetting analysis (10) confirm KOI-872.01 and the perturber 
as being genuine planets, henceforth referred to KOI-872b (corresponding to KOI-872.01) and 
KOI-872c. 

The period ratio, P c / 'P b = 1.697 indicates that the two planets are just outside the 5:3 reso- 
nance. This may be a relatively common configuration probably related to the radial migration 
of planets in the protoplanetary nebula (23). The resonant angle, 5A C — 3A& — 2w c , circulates in 
the retrograde sense with the period of ~ 2 yr. The dominant TTV period (~ 5.6 transit cycles; 
Fig. |2]) comes from the relatively distant 2: 1 and 3:2 terms (periods ~ 190 yr). 

The orbital eccentricities are e& < 0.02 and e c ~ 0.015. The nearly circular orbits probably 
indicate that the two planets formed at, or migrated to, their present orbits without suffering 
any strong dynamical instability. The two planets also exhibit nearly but not exactly coplanar 
orbits. The orbital inclination of KOI-872c with respect to the transit plane is J c < 4.5° (99% 
confidence interval) with the best TTV fit indicating I c = 2.6°. The best-fit inclination value 
and fi c ~ 298° obtained from the TTVs suggest the planet just avoids transiting. A search for 
the transits of KOI-872c indeed yielded no events, implying that I c > 1°. 

Although no transits of KOI-872c were detected, we do detect a 9-cr transit signal on a short 
period of 6.8 d. The transit corresponds to a 1.7 R® Super-Earth, which we refer to as KOI- 
872.03. Because of the probable low-mass nature of this body, the expected TTVs induced on 
KOI-872b would be around 1 s in amplitude, too small to detect with our data. Further, the 
TTVs of KOI-872.03 itself are estimated to be ~ 10 s, also undetectable. Without TTVs, we are 
unable to confirm that KOI-872.03 orbits the same star as KOI-872b and KOI-872c, as opposed 
to a blended background star. However, our analysis finds that KOI-872.03 's light curve derived 
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stellar density is consistent with that of KOI-872b for both planets on near-circular orbits. 

We predict that the radial velocity (RV) measurements of KOI-872 should reveal at least 
two basic periods. The RV term with a 57.0 d period, corresponding to KOI-872c, should have 
K c ~ 20 m s _1 half- amplitude. The amplitude of the 33.6 d period term, corresponding to 
KOI-872b, is uncertain because we do not have a good constraint on M b . The 6.8 d term, 
corresponding to KOT872.03, will be difficult to detect because K 3 ~ 1-2 m s^ 1 (assuming 
Earth-like density). 

To gain insights into the system's dynamical behavior we numerically integrated the plane- 
tary orbits starting from the best-fit solution (Fig. [3]). The semi-major axis of KOT872b expe- 
riences short-period variations, related to the TTVs, with an amplitude of ~ 5 x 10~ 4 AU. The 
eccentricities undergo anti-correlated oscillations, as dictated by the angular momentum con- 
servation, around means e b = 0.013 and e c = 0.012, with a period of ~200 yr (corresponding to 
zu b s precession period; zu c precesses much slower). The anti-correlated oscillations of inclina- 
tions relative to the transit plane are a geometrical effect resulting from the shared precession of 
£l b and fi c (~ 200 yr period). The inclinations relative to the invariant plane are nearly constant 
with the mutual inclination ~ 1°. If I b relative to the transit plane increases in the next decades 
(Fig. 3d), KOI-872b's transits will gradually disappear. 

TTVs were originally proposed as a nontransiting planet detection method (2-4), but have 
recently found more use in validating the transiting planet candidates from Kepler (24, 25). 
Kepler has previously inferred the presence of a nontransiting planet via TTVs, but showed that 
the measurements were unable to support a unique solution (26). Here we have demonstrated 
the full potential of TTVs as a method to detect nontransiting planets and precisely characterize 
their properties. 
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Table 1: KOI-872 system parameters. KOI-872b parameters were computed from the weighted 
posteriors of a model accounting for TTVs, using MultiNest. Parameters fitted in the transit 
model are quoted as the median of the marginalized posteriors with ±34.13% credible intervals. 
Instrumental terms and times of transit minimum may be found in Table S3. KOI-872c param- 
eters were computed from the fit to the KOI-872b's TTVs. Parameters fitted in the TTV model 
are quoted as maximum likelihood with ±34.13% uncertainties computed using the A% 2 = 1 
method described in (27), where the TTV errors have been rescaled such that Xreduccd = 1- 
The 99% confidence areas from the TTV fit alone are shown in Fig. S10. The measured TDVs 
do not offer a meaningful constraint, since the parameter sets near si that fit the TTVs also 
fit the TDVs. Orbital longitudes of KOI-872c are relative to the transit reference system on 
2455053.2839 BJDutc- KOI-872.03 parameters were computed from an MCMC run. Moon 
mass constraint (3 a limit) was derived from model Mmt2,ro {10). 
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Figure 1 : Maximum likelihood transit model (red line) overlaid with the long-cadence Kepler 
offsetted data for KOI-872b. The large TTVs are evident visually from the light curve. The 
ramp-affected transit is excluded here (see Fig. S2). 
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Transit Cycle 




Transit Cycle 



Figure 2: Transit timings variations. The measured TTVs (from model M.t) and TDVs (from 
model Aiy) and their uncertainties are indicated in panels (a) and (c). Panel (a) shows the 
calculated values for si (red line). The TDVs of si are consistent with the measured, flat TDV 
profile. Panel (c) shows s2 (blue line) and our best moon model (green line). While s2 also 
fits the measured TTVs relatively well, the strong TDV trend in (c) is inconsistent with the 
measurements. Panels (b) and (d) show our predictions for si and s2. 
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Figure 3: Orbit evolution. Figure shows the semi-major axis (panels a and c), eccentricity (b) 
and inclination (d) for the best fit solution. The orbital elements of KOI-872b and KOI-872c are 
shown in red and blue, respectively. The inclination in (d) is defined relative to the transit plane 
(I = 0° would correspond to a transit across the middle of the host star's disk). An approximate 
transit zone for Q = 270° is shaded. and fi c are locked near 270°, because the invariant plane 
of planets is tilted to the transit plane. KOI-872.03 is omitted here because it has a negligible 
effect on the dynamical evolution of KOI-872b and KOI-872c. 
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Supporting Online Material 



System Identification We here discuss why KOI-872 was selected for analysis by the HEK 
project. The HEK project filters the list of all known Kepler transiting planet candidates down 
to a subset of the most promising for detecting exomoons. This process is described as "Tar- 
get Selection" (TS) and features three distinct pathways: target selection visual (TSV), target 
selection automatic (TSA) and target selection opportunities (TSO). We direct the reader to (8) 
for details on each procedure. 

KOI-872b, internally labeled as HCV/HCA-439.01, was identified independently by TSV 
and TSA as being an excellent candidate for exomoon follow-up and thus was prioritized for 
more detailed analysis. These determinations were initially based upon the Q1-Q3 long-cadence 
data only. Subsequent Q4-Q6 data is included in the analysis presented here. 

Data Handling We will here describe the sequential steps we took in processing the Kepler 
photometry for the target KOI-872. 

Data acquisition. We make use of the publicly available archival data from the Kepler 
Mission on the MAST website, which consists of quarters 1 through 6 (Q1-Q6). All data is 
in long-cadence mode with nearly continuous coverage. Full details on the data processing 
pipeline can be found in the data release handbooks and is not repeated here. 

Long-term detrending with a cosine filter. We make use of the "raw" (labeled as "SAP_FLUX" 
in the header) data processed by the Kepler DAWG (Data Analysis Working Group) pipeline. 
All of the publicly available photometry was acquired in long-cadence (LC) mode spanning 
quarters 1 through 6 (Q1-Q6) and data releases 4 through 9 (DR4-DR9). A detailed descrip- 
tion of the pipeline can be found in the accompanying release notes. The "raw" data has been 
processed using PA (Photometric Analysis), which includes cleaning of cosmic ray hits, Ar- 
gabrightenings, removal of background flux, aperture photometry and computation of flux- 
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weighted centroid positions. 

The data releases also include corrected fluxes (labeled as "PDCSAP_FLUX" in the header), 
which are outputted from the PDC (Pre-search Data Conditioning) algorithm developed by the 
DAWG . As discussed in DR5, this data is not recommended for scientific use, owing to, in part, 
the potential for under/over-fitting of the systematic effects. For the sake of brevity, we do not 
reproduce the details of the PA and PDC steps here, but direct those interested to (28) and the 
DR handbooks. 

The Ql to Q6 PA photometry are shown in Fig . [T(a)| [T(b)l [T(cjl [T(d)| [T(e)1 &[Tffl1respectively. 
One challenge in attempting a correction is assessing which components are astrophysical in 
nature and which are instrumental. The astrophysical signal of interest is the transit signal and 
thus we aim to detrend all other effects and preserve the eclipse. The spurious trends may be 
removed by applying a high-pass filter to the photometry, in a similar way as was used by (29) 
for CoRoT photometry and (30, 31) for Kepler photometry. To remove the long-term trends, we 
thus applied a discrete cosine transform (32) adopted to the unevenly spaced data. 

We first removed the transit events with a margin of ±Xdays either side of the times of 
transit minimum, as predicted from the linear ephemeris period and epoch reported in (9). The 
choice of the time T is provided later. We also remove outliers, identified as those points lying 
3-cr away from a spline-interpolated running median of window-size 600 minutes. Treating 
each quarter separately, we fitted the remaining data with a linear combination of the first N 
low-frequency cosine functions using 



where tj is the timing of the j measurement, i = 0, iV in integer steps and N is equal to 
the rounded integer value of (25/41') where B is the timespan of the observations and X' is 
the timescale we are protecting. Since we are trying to protect the timescale X, we set V = 




(1) 
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3T to ensure the timescale of interest in minimally distorted^. Using a Levenberg-Marquardt 
algorithm, we then fit for the linear coefficient, for each of the cosine functions, so that the 
fitted model is given by 



N 



2R(*i) = I>/i(*i)- (2) 

i=0 

We then subtracted model 9Jt from the light curve (including the transits). The model is 
shown over the data for the Q1-Q6 photometry in Fig. [1(a)} [T(b)l [T(c)l [T(d)l [T(e)1 & [T(f)} The 
final time series is clipped to be ±21 surrounding the time of transit minimum, as based upon 
the linear ephemeris of (9). 

Protected timescale. The harmonic filter protects a chosen timescale. When searching for 
ellipsoidal variations or reflected light, this timescale may be set to be equal to the period of the 
transiting planet, as was done in (29) and (33). However, if we are solely interested in the transit 
event, then the timescale of interest is much shorter than this. Specifically, it is the duration of 
the transit which is the timescale of interest. If we are looking for moons, then auxiliary eclipse 
signals may appear either side of the transit event. Thus, we wish to protect a timescale which 
includes both the transit duration and the surrounding temporal coverage corresponding to the 
region where a moon could reside. Spatially, this is the Hill radius and thus we need to define 
"a Hill timescale". 

The first-to-fourth transit duration is already known from (9) to be T i4 = 4.3863 hours. If 
this time is roughly equal to the time is takes for the planet to traverse two stellar radii, 2/?*, 
then the velocity of the planet is given by vp ~ 2R*/Ti4. The Hill timescale, T#, is given by: 



'For Q6, we were not able to get a good fit to the raw data using T' = 31 and so this factor was relaxed to 
unity. 
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apTu I Mp 



1/3 



2R* \3M* 



(3) 



where R H is the Hill radius, a P is the planet's semi-major axis, M P is the mass of the planet 
and M* is the mass of the star. In the above, a P , T i4 , R* and M* have reasonable estimates 
from (9). M P , however, is wholely unknown. Since the planet is the Saturn-size regime, we 
simply adopt a Saturn-like density and estimate M P accordingly: 



where Rp and pp are the planet's radius and mean density respectively. For HCV/HCA- 
439.01, we estimated T H ~ 0.24 days. The protected timescale, T, was then set to be equal 
to: 



where the 1.2 factor is added as a 20% buffer. In total, this gives X = 0.80 days as the 
protected timescale. 

Dilution factors. In each aperture, a small amount of third light is typically present from the 
overlapping point- spread-functions in Kepler's relatively crowded field. To account for this, it 
is necessary to include the third light in the transit light curve fits. We follow the method of (34) 
to make this correction. Although not available in the MAST headers directly, the dilution 
factors are accounted for in the PDC photometry but not in the PA photometry. Therefore, we 
simply take the median of both sets to compute the dilution factor. This is done for each quarter 
independently and ranged from 3% to 9%. 




(4) 



Z=1.2(T U + 2T H ), 



(5) 
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Light Curve Analysis Star spots. The raw transit light curve shows clear evidence for rotat- 
ing star spots on the surface of KOI-872. In particular, Fig. |l(d)| shows three deep minima due to 
a large dark, rotating spot of roughly 10% the radius of the star. However, throughout the time 
series complex spot patterns emerge and evolve indicating the behavior is more complex than 
a simple, single spot model. For this reason, a simple periodogram is unable to reliably infer 
a rotation period. A more complex spot-model, accounting for multiple spots and differential 
rotation, is required to interpret the behavior of the these events, which is outside of the scope 
of this work. 

For the purposes of this work, the most important consequence of the spots is that signals 
which appear like exomoon mutual-events may in fact be star spot crossings. Equipped with 
this prior information, a star spot crossing seems an a-priori more likely explanation for the 
TSV signals initially identified. These signals also complicate the derivation of upper limits 
on a putative exomoon, as discussed later. Upon downloading the Q4 data in early January, 
we immediately suspected the putative signal was a false-positive. However, we continued our 
analysis because our early model fits revealed that the planet exhibited large and complex transit 
timing variations, as also reported in (35, 36). 

Model fits. In fitting the transit light curve, we consider a variety of models to explain 
the data. In all cases, the fitted parameter set is the same as that described in (8) with two 
exceptions. Firstly, rather than fitting for the satellite-to-planet mass ratio (Ms/M P ) in planet- 
with-moon fits, we used (ps) 2 ^ 3 i-e. mean density of the satellite to the power of two-thirds. 
This was done since a physical density has a better known prior than mass ratios. The two-thirds 
power ensures uniform priors in any derived mass-ratios since the other density terms feature 
this index for reasons described in (8). Secondly, we fit for a photometric noise term, aw, 
which is the standard deviation of the noise assuming it to be white. This is done to propagate 
the uncertainty of the photometric uncertainties themselves into our calculation of the Bayesian 
evidences, aw is fitted with a modified Jeffrey's prior with the inflection point at the median 
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photometric uncertainty from the DAWG pipeline and the maximum limit being 10 times larger. 
We note the best value for this term (from the Mt fits) is aw = 2.603^o;o2i mmag per minute. 

In what follows, we assume that the system is a genuine planetary system rather than a 
false-positive, such as a blended eclipsing binary. Details on our vetting procedure and blend 
analysis are provided later. The following principal models were considered: 

■ Mp - Planet- only 

■ M.t - Planet-only with variable times of transit minimum (TTV) 

u M.y - Planet-only with variable times of transit minimum (TTV), transit depths (T<5V), 
impact parameters and [p" rc ] 2 ^ 3 values^ (thus permitting TDV) 

■ M.m - Planet- with-moon 

■ M.m,rq - Planet-with-moon, defining the moon as a point-mass i.e. we fix Rs = 0. 

■ -M. mti,mo ' Planet-with-moon, removing the maximum a-posteriori transit timing varia- 
tions deduced from Air- Since all TTVs are removed, we must enforce Ms = 0. 

■ M. MT2,m - Planet-with-moon, removing the maximum likelihood transit timing varia- 
tions from a second planet fit and searching for residual TTVs/TDVs only. 

The fits were executed using the MultiNest algorithm (11, 12), which is a multimodal 
nested sampling routine (37) designed to compute the Bayesian evidence in complex parameter 
space in an efficient manner. For brevity, we direct those interested to the aforementioned works 
for further details. MultiNest is coupled with the forward-modeling code of LUNA (38), 
which is designed to model the transit light curves of a planet-with-moon accounting for mutual 
events, auxiliary transits, dynamical perturbations and non-linear limb darkening in an analytic 
manner. 

2 p" IC is the light curve derived stellar density assuming a circular orbit. 
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Quadratic limb darkening coefficients were estimated in a similar manner to that described 
in (30). For this calculation, we assumed the effective temperature and surface gravity of the 
star to be that reported in (9) (T eff = 5127 K and log g* = 4.59), which in turn come from 
the Kepler Input Catalogue (KIC). Our later spectroscopic analysis shows these to be excellent 
estimates. For the Kepler bandpass, we used the high resolution Kepler transmission function 



found at http://keplergo.arc.nasa.gov/CalibrationResponse.shtml We employed the atmosphere 
model database from (39) providing intensities at 17 emergent angles, which we interpolated 
linearly at the adopted T c g- and log g* values. The passband-convolved intensities at each of 
the emergent angles were calculated following the procedure in (40). This whole process is 
performed by a Fortran code written by I. Ribas. To compute the coefficients we used the limb 
darkening law given in Equation [6] 



Jt = i _ _ ^ _ U2 (i - ^ (6) 
H 

where the various terms are defined in (40). The final coefficients resulted from a least 
squares singular value decomposition fit to 1 1 of the 17 available emergent angles. The reason 
to eliminate 6 of the angles is avoiding excessive weight on the stellar limb by using a uniform 
sampling (10 fi values from 0.1 to 1, plus ji = 0.05), as suggested by (41). This leaves us with 
Mi = 0.3542 and u 2 = 0.3607. 

Ramp correction. During the Kepler time series, there are a few safe mode events where 
the telescope stopped observing. These result in a pause in the continuous photometry followed 
by an exponential ramp as the telescope starts observing again. This ramp, likely similar to the 
charge trapping effect seen with Spitzer (42), lasts for around one week and we usually simply 
clip the affected data. Unfortunately, the second transit observed occurs during one these ramps. 
In order to maximize the available data, we decided to correct for the ramp effect and recover 
this transit. We apply a simple exponential decay model of the form 
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F'/F = a - a-y exp(-t/a 2 ), (7) 

where F' is the flux affected by the ramp, F is the flux corrected for the ramp and t is the 
time since the start of the ramp. We also tried a double exponential, similar that advocated 
by (42), but found the two timescales converged to a single value. In fitting the data, we fit the 
parameters a (absorbed by the OOT vector in practice), a\ and a 2 simultaneously to fitting 
the transit model. We find that this simple model provides an excellent description of the ramp 
effect, as visible in Fig [S2j We note that the ramp timescale was found to be best modeled by 
d2 = 1. 10241^0069 days (from _Mt fits), which may be useful to other observers. 

Model selection. In searching for an exomoon, one must conduct model selection between 
the various hypotheses which could explain the data. We define our null model to be that of a 
transiting planet without a moon and with static parameters; model M. p. This simply assumes 
a constant linear ephemeris with a constant duration and depth every transit. 

We also consider models of a planet without a moon, but with perturbations. The simplest 
type of perturbation model we consider is that of a planet with varying times of transit minimum 
(i.e. a planet experiencing TTV) which we dub M. T . Model M. v extends this to allow for 
variable depth and duration as well. Models M.m and A4m,ro are the planet-with-moon model, 
as simulated from LUNA, except the latter assumes a fixed zero-radius moon (R s = 0); i.e. 
considers TTVs and TDVs only. 

Due to the presence of a second planet inducing the TTVs (as discussed in the main text), we 
also tried models Ai mti,mo and M.mt2,ro- The first model removes all transit timing variations 
by subtracting the maximum a-posteriori transit times deduced using M. T from the original data 
and then refitting for a zero-mass exomoon (we cannot fit for a moon mass if we have forcibly 
removed all TTVs). The second subtracts the maximum likelihood model for a second planet 
(discussed later) causing the TTVs and then fits the adjusted data for a planet-with-moon. In 
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this case, we also fix the moon radius to be zero since we later hypothesize that stellar activity 
may be inducing false-positive moon-like eclipses. Therefore, this model looks for residual 
TTVs and TDVs only. Because these fits use modified data for inputs, they cannot be directly 
compared to the other models and require a custom null-model for comparison in each case. 

Model selection is performed by comparing the Bayesian evidence of each model. The 
higher the Bayesian evidence the more likely the model is the correct one. Computing the 
Bayesian evidence is computationally expensive and particularly challenging in the high-dimensional 
space we are faced with. For example, Mr involves 38 free parameters. 

In some cases, we found it was not necessary or possible to formally compute the absolute 
value of the Bayesian evidence, Z. For example, as MultiNest iteratively searches through 
higher log Z values, the code can be stopped if the log Z value greatly exceeds the evidence of a 
competing hypotheses. This allows us to place a lower limit on the confidence of such a model, 
as was done for M. T . Here, we found (\ogZ T — \ogZ v ) > (166.1 ± 0.7) (where model My 
has the next-best Bayesian evidence), indicating > 18.1-er preference for the TTV model over 
the TTV+TDV+T5V model. To highlight the computational demands, we point out that even 
this lower limit required over 3 years of equivalent processing time, including over 4 billion 
likelihood evaluations, with a 2. 1 GHz AMD Interlagos CPU. We also note that the TTV model 
is preferred over the static model, Mp, at a confidence of > 43.9-cr, which to our knowledge 
represents the highest formal significance for a TTV detection ever reported. For cases where 
only a lower limit on Z is provided, the posteriors were computed by re-running MultiNest 
in constant efficiency mode. 

Models Mm and Mm,ro are found to be relatively poor fits to the data with (log^M — 
\ogZ T ) < -(635.3 ± 0.4) and (logZ MjR0 - \ogZ T ) < -(545.8 ± 0.4). This is visually 
evident by comparing the TTV model for a second planet versus that of a moon in Fig. 1 . Since 
Z T 3> Z P , there is a very high probability that transit timing variations are present. But since 
a moon model provides a much lower evidence than the TTV model (i.e. Z M «C Z T ), then we 
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deduce that i) TTVs are present ii) a moon is not responsible. 

Although this conclusion tells us the large TTVs are not being caused by a moon, they do 
not exclude the presence of a moon either. The possibility of a moon as an independent source 
of TTVs is investigated using models Mmti,mo and Mmt2,ro- We did not directly compare 
the Bayesian evidence for these models with the others. This is because the input data was 
manipulated in each case by subtracting times of transit minimum from each epoch. In order to 
perform a reliable model comparison, we generated a custom null model for each: M. Mri.nuii 
and M.MT2,xm\\- These were accomplished by re-running the fit with the size and mass of the 
moon set to zero, allowing us to remove the moon terms as free parameters. The results of these 
fits will be discussed later in a dedicated section. 

The priors and Bayesian evidence values for each model we attempted are provided in Ta- 
bles [STT&fS2l respectively. Table 1 provides the final physical system parameters. Table [S3] 
provides the final system parameters for instrumental terms and times of transit minimum. 

A Search for an Occultation Searching for an occultation is challenging due to the large 
transit timing variations present. In order to accommodate for this, we allow the occultations 
to have their own timing variations. To minimize the number of free parameters, the p, bp, 
[p* irc ] 2 / 3 , P P and r parameters were sampled from a Gaussian prior derived from the M. T 
model fit. Since any eccentricity induced timing offsets should be absorbed by the occultation 
timing variations, fitting for eccentricity terms would be essentially fitting redundant parameters 
and thus we fix ep = 0. This left us with 15 OOT parameters, 15 occultation times, 1 noise 
parameter (aw) and 1 term for (J^/J 7 *), the flux-per-unit-area ratio of the planet and star. We 
ran two fits; one where (J-p/J 7 *) = (in such a case the occultation times are not required 
as free parameters) and one where < (J-p/J 7 *) < 1 and was a uniform prior, allowing us 
to perform a model comparison later. Note that under the assumption that the system is a real 
planetary system, the TTVs suggest ep ~ anyway, as described in the main text. 
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A comparison of the Bayesian evidence from the null fit versus the occultation fit yields 
(log Z occ — log Z im n) = (—2.43 ± 0.44). Thus, the null hypothesis of there being no occultation 
present is the preferred model. 

Although no occultation is detected, we may use our results to place upper limits on the oc- 
cultation depth. The marginalized posterior of of the occultation depth yields 5 occ = 9.9+7 4 4 9 ppm, 
and places a 3-a upper limit of S occ < 71.0 ppm. Assuming a geometric albedo of unity 
(A g = 1), the reflected light component of the occultation is expected to be 3.9 ppm. As- 
suming the occultation is due to reflected light only, our occultation depth limit corresponds to 
A g < 18.0 i.e. we are unable to constrain the albedo of KOI-872.01 to any physically plausible 
range. 

Given our insensitivity to reflected light, our upper limit is more robust and meaningful 
when interpreted as an upper limit on thermal emission. Kepler's visible bandpass is not well- 
suited to detecting thermal emission but a meaningful constraint can still be derived. Treating 
the star and planet as black bodies and integrating over the custom Kepler bandpass, we find 
T P < 2442 K to 3-a confidence (assuming T* = 5155 K). 

Vetting the Planetary System So far, we have assumed the observations are due to an un- 
blended planet transiting a star. Here, we consider alternative models which do not require the 
system to include a planet. To avoid confusion with the earlier model fits which assume the 
eclipsing object is an unblended planet, we will dub these models as hypotheses, Hi. There at 
least three such hypotheses which could potentially explain the data: 

■ "Hp: No blend is present and thus we have a planet transiting a star 

■ Heb, 33.6 : A blended eclipsing binary (EB) with the eclipsing bodies on an orbital period 
of 33.6 d 

■ Heb,67.2'- A blended eclipsing binary (EB) with the eclipsing bodies on an orbital period 
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of 67.2 d 



In the last two cases, the eclipsing binary could be a larger planet (e.g. Jupiter-sized) eclips- 
ing a star and thus could still be considered a genuine planetary system. Further, the blend 
source could be foreground, background, associated or a mixture of multiple sources. The last 
two cases can also be considered in two flavors i) EB on a circular orbit ii) EB on an eccentric 
orbit. The possibility of an unblended grazing eclipsing binary is included in the model Hp, 
and a blended grazing eclipsing binary within V,eb,33.g and Heb, 67.2- 

In general, one expects the false-positive rate for Kepler Objects of Interest to be quite low, 
with recent estimates arriving at < 10% (43). Nevertheless, we will here investigate the possi- 
bility of the blended EB scenarios mimicking a planetary system. We will approach this problem 
using several tools 1) a centroid analysis 2) constraints from the spectroscopy 3) model selec- 
tion with Bayesian evidence determinations of the transit light curve shape (a blend analysis) 4) 
dynamical constraints from the timing variations and stability arguments. 

A Centroid Analysis Overview. The DAWG pipeline output provides flux- weighted centroid 
positions for all observed targets. (44) have demonstrated that the very small shifts in centroid, 
expected for blended occulting sources, can be measured accurately from the Kepler data. Con- 
sider two physically separated sources with overlapping PSFs. The computed centroid position 
is flux- weighted with respect to these two sources. When one of the sources is eclipsed, its 
flux temporally decreases and thus the flux-weighted centroid shifts towards the other source. 
Therefore, the detection of a shift in flux- weighted centroids during the eclipses would indicate 
the presence of a blend source. An example of this technique detecting such a source is for 
KOI- 13 (45). 

Excluded Centroid Shift. We extracted the x and y centroid positions for KOI-872 surround- 
ing ±0.4 d of the predicted transit events according to the linear ephemeris model derived in 
M p. We then removed the maximum a-posteriori times of transit minimum, r, for each transit 
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epoch, as computed by the M. T model. This step essentially phases all of the data. Finally, we 
then divided epoch by the median x and y centroid position to remove the effect of long-term 
trends in the centroid positions. No cleaning or detrending of the centroids was attempted. 

The phased centroid positions, shown in Fig. [S3l display no obvious up or down pattern 
during the transit events (marked by the vertical gridlines). The scatter can be seen to be com- 
parable, but somewhat larger than, that derived for KOI- 13 in Fig. 2 of (45). Taking the mean 
and standard deviation of the in- versus out-of-transit centroid positions we find: 



Ax = x out - x in = (0.09 ± 0.57) x 10" 4 , (8) 
Ay = y out - y m = -(0.07 ± 0.58) x 10" 4 . (9) 

Defining Ar = a/ Ax 2 + Ay 2 , we determine Ar = (0.12 ± 0.57) x 10~ 4 , where all units 
thus far have been given in units of pixel position. This analysis clearly indicates that the 
data are consistent with no separated blend source present. The 3-er upper limit, converted 
to arcseconds, corresponds to Ar < 0.68 mas, demonstrating the impressive performance of 
Kepler once again. 

Excluded Blends. We here describe a toy model to interpret this upper limit. We consider 
two sources of flux and F B (star and blend source), where the star is transited. The flux- 
weighted positions, r, of the centroid in- and out-of-transit are given by: 



r*F* + r B F B ^ 



F* + F B 
r*F*(l - 6) + r B F B 



(11) 



F,(l-5) + F B ' 

where 5 is the unblended eclipse depth. Since we only detect the blended eclipse depth, <5 b s , 
we must convert between the two using the expression from (34): 
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5 



<Wl + /3), 



(12) 



where f3 = F B /F*. If we also make the replacement r B = + Ar, one can write: 



^'"obs '"in '"out > 

(35 ohs Ar 



(13) 



(14) 



Here, Ar obs represents the observed change in centroid position during the eclipse, in arc- 
seconds, and Ar represents the physical separation of the two sources on the sky, in arcseconds. 
Defining the blend factor as B — (1 + 0), one may solve inverse the above expression to solve 
for B, as a function of Ar: 



Using this simple model, we plot the excluded values of B as a function of Ar in Fig.lS4l 
Note that this result is purely based on the centroid shifts. To quote several values from the 
figure, B < 1.239 for Ar < 0.5", B < 1.051 for Ar < 2" and B < 1.016 for Ar < 6". A 
blend factor of ~10% does not impact significantly on our results and thus only blends within 
Ar < 0.5" could possibly cause a planet false-positive. We note that, in general, such a closely- 
space companion is quite rare, with a recent adaptive optics campaign on suitable KOIs finding 
only 6.7% of KOIs have a companion within 0.5". 

Fig. |S4] also provides the same constraints as computed for the inner transiting planet candi- 
date KOI-872.03 (dashed line). Details on the detection of this candidate are provided later. Due 
to the much smaller depth, the upper limits are less constraining but nevertheless still consistent 
with the absence of a blend source. 



B(Ar) 



<5ob s Ar 



(15) 



<5ob s Ar - Ar obs + <5 obs Ar obs 
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A Spectroscopic Analysis. Observations. Spectroscopic observations of the star HCV/HCA- 
439 were carried out with the Astrophysical Research Consortium Echelle Spectrograph (ARCES) 
on the Apache Point Observatory 3.5 m telescope located at Apache Point Observatory (APO) 
in New Mexico. We used a 176 x 372 slit which delivers a spectral resolution of AA /A 3 1 ,000 
over a spectral range of 3200-10,000 A. We obtained a total of two thirty minutes exposures on 
the star HCV/HCA-439 on the night of 2012 Jan 31 which, when combined, yielded a total S/N 
of 11 at 5100 A. After performing a standard overscan correction, we removed cosmic rays, 
extracted the spectra, applied a flat-field correction, and determined the dispersion correction 
from a ThAr lamp spectrum using standard techniques in the IRAF's IMRED, CRUTIL, and 
ECHELLE packages. 

Stellar Parameters. We used the Yonsei-Yale (YY) isochrones (46) to determine the phys- 
ical properties of the host star. The first input into this analysis were the stellar atmosphere 
parameters determined from the APO 3.5 m spectra. We used the Stellar Parameter Classi- 
fication (SPC) method to derive the stellar atmosphere parameters. SPC cross-correlates the 
observed spectrum against a grid of synthetic spectra drawn from a library calculated by John 
Laird using Kurucz models (47). The synthetic spectra cover a window of 300 A centered near 
the gravity-sensitive Mgb features and has a spacing of 250 K in effective temperature, 0.5 dex 
in gravity, 0.5 dex in metallicity and 1 km s _1 in rotational velocity. To derive the precise stellar 
parameters between the grid points, the normalized cross -correlation peaks were fitted with a 
three dimensional polynomial as a function of effective temperature, surface gravity and metal- 
licity. This procedure was carried out for different rotational velocities and the final stellar 
parameters were determined by a weighted mean of the values from the spectral orders covered 
by the library. 

The isochrone analysis made used of the stellar effective temperature T cff = (5 155 ± 105) K, 
and the metallicity [Fe/H] = (0.41 ±0.10) from the SPC analysis. The second input into the YY 
analysis came from the light curve analysis of the Kepler data. The transit duration is closely 
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related to the ap/R* parameter, which in turn determines the mean stellar density p* (48). In 
general, this trick is only possible for systems where the orbital eccentricity is precisely known 
(49). It will be shown in the dynamical analysis discussion later that planet b's eccentricity is 
indeed strongly constrained from the TTVs to be near-circular. We proceed by using the light 
curve derived stellar density posterior from our model fits, giving = 1520^^0 kg m~ 3 . 

The mean stellar density acts as a luminosity indicator for the star; smaller density typically 
means more evolved stars. Another possible luminosity indicator would be the surface gravity 
of the star, as determined from the spectroscopic analysis. If the eccentricity of the system is 
well determined and the light curve is of high quality (both of which are true here), then ap/R* 
and the corresponding is typically a better luminosity indicator than log g*, in the sense that 
the derived physical parameters have a smaller error. We have generated over 10,000 values of 
T c g, [Fc/H] and ap/R* using their a posteriori distribution (assuming Gaussian distribution for 
the first two), and searched the YY isochrones for each materialization. About 95% of the input 
values had a matching isochrone. Stellar parameters were then determined as the median of the 
resulting distribution. The final parameters were M* = 0.90 ± 0.04 M , R* = 0.94 ± 0.04 R Q , 
and logg* = 4.44 ± 0.04 (cgs). The isochrones and the final solution are shown in Fig. IS5l 
where the backdrop of isochrones is for ages 0.2, 0.5, 1 .0, 2.0, . . . 13.0 Gyr (from bottom to top) 
and [Fe/H] = 0.41. The solution indicates an old star with 10 ± 3 Gyr age. 

In regard to vetting of the system, we find no evidence for double-lines indicating a blend. 
We also find the spectral classification of the star very well described as a dwarf main-sequence 
star, and exclude the possibilities such as a giant star or white dwarf. We note that our classifi- 
cation is in close agreement with the KIC determination reported in (9). 

A Blend Analysis Overview. The third vetting tool we use is model selection with Bayesian 
evidence determinations of the transit light curve shape i.e. a blend analysis. Blend anal- 
yses have become a powerful instrument in the toolbox of the Kepler team too, using their 
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custom BLENDER software (50, 51). Here, the team simulate a grid of billions of possible 
false-positive scenarios and compute the odds ratio of valid planet solutions versus valid false- 
positive solutions. Systems strongly favoring the planet solution are considered to be "vali- 
dated". BLENDER makes use of additional information such as centroid positioning, priors of 
the frequency of eclipsing binaries and multi-color light curves (often from Spitzer (51)). In 
this work, we limit our analysis to purely an inspection of the shape of the light curve. This is 
done to reduce the computational demands of simulating billions of false-positives, and yet take 
advantage of the fact we have other more powerful constraints from the dynamics (as discussed 
later). Further, a centroid analysis has already indicated that a blend must be within 0.5" to have 
a significant impact on our results and thus can be treated as a separate line of evidence. 

Although our blend analysis is more simplified than BLENDER, it does take advantage of 
the Bayesian evidence for model selection, which is currently not implemented in BLENDER 
or other blend analyses in the current literature. Our approach is to consider a null hypothesis 
and then more elaborate models involving blend scenarios and compute the Bayesian evidence 
of each model. Since the more elaborate models include more parameters (i.e. a greater prior 
volume), they are penalized in the computation of the Bayesian evidence. If a blend model 
provides a Bayesian evidence significantly greater than that of the null model (A \ogZ > 5), 
it passes the first test to becoming the preferred model. The second test we impose is that the 
parameter posteriors of the model must correspond to a physically plausible scenario. If both 
of these criteria are satisfied, then the null model can be displaced. 

Another possibility is that two or more hypotheses yield approximately equal Bayesian evi- 
dences, such that there is no statistically significant preference between them. In such a case, we 
continue to consider these hypotheses as plausible descriptions of the system and test whether 
they are consistent with our subsequent dynamical analysis too. 

The null hypothesis is the simplest model which can explain the data and so involves the 
fewest parameters to describe the system. In our case, this represents just two eclipsing objects 
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without a blend. Whilst an unblended grazing eclipsing binary could fall into this category, it 
will be shown shortly that such a hypothesis is highly improbable. Therefore, the null hypothe- 
sis essentially represents a planet-sized object transiting a star. 

Implementation. In our blend analysis, we are only investigating the shape of the light 
curve. The dynamical constraints from the timing variations will be discussed later, and so 
here we eliminate them by subtracting any TTVs away from a linear ephemeris. This is ac- 
complished using the maximum a-posteriori transit times from the model fit M. T performed 
earlier. This reduces the number of free parameters by 15 and makes the fits far easier to handle 
computationally. However, we still fit each eclipse epoch with a unique baseline to remove any 
residual DC power from the detrending procedure. 

For the null hypothesis, "Hp, we have 15 OOT baseline parameters, p, bp, [p* irc ] 2 / 3 , Pp, r , 2 
instrumental terms to describe the ramp effect for the one affected event and two limb darkening 
parameters (u\ and [u% + w 2 )), giving 24 parameters in total. Limb darkening was fitted for to 
provide a fair comparison to the blend models where limb darkening cannot be assumed to be 
the same as the theoretical models used in the planetary fits. The priors on the limb darkening 
terms were selected such that the brightness profile is positive everywhere and monotonically 
decreasing from limb to center, specifically < U\ < 2 and < (u\ + u 2 ) < 1 (52). 

We also extended the prior on p to the range < p < 1 and similarly for the impact 
parameter < bp < 2. The other fitted terms had the same priors as used before (see Table [STT). 

For the 33.6 d period EB scenarios, we simply add a single additional term, a blending 
factor B. Our model for the diluted light curve follows the prescription of (34), where B = 1 
indicates no blend and B > 1 indicates a blend. The model is general in that the source of 
the blend could be foreground, background, associated or a mixture. The prior on the blending 
factor was chosen to be 1 < B < 100, and all other priors were left unchanged from "Hp. In 
general, one expects the 33.6 d blended EB scenario to give an equally good fit to the data as a 
planet, since the B term is degenerate with the other fitting parameters (34). For the eccentric 
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cases, we fitted for e P and u P directly using < u P < 2ir and < e P < 0.9 uniform priors. 
The eccentricity was cut-off from extremely eccentric orbits to save CPU time, since solving 
Kepler's equation takes dramatically longer in the extreme eccentricity regime. 

The final scenario of a 67.2 d period blended EB was treated by first modifying the prior 
on the period to be uniform around ±1 day of 67.2 d. Then, we instructed the code to treat the 
occultations has having zero-limb darkening. Whilst the ingress/egress of the occultation may 
in fact have limb darkening, the 2nd-to-3rd contact of the transit, which dominates the signal- 
to-noise, cannot have a limb darkened profile. Therefore this approximation contains the most 
important physics of the problem. The primary transit is treated as a limb darkened event as 
before. The occultation depth is equal to p 2 (J r P / J 1 '^), where T is the flux-per-unit-area of a 
body. This flux ratio is the only new parameter required, for which we use < (J^/J 7 *) < 2 
as our prior. In general, one expects this scenario to be easier to distinguish against a planet due 
to the lack of curvature in alternate eclipses. 

Results: An Unblended Grazing EB. The easiest scenario to disregard is that of an unblended 
grazing eclipsing binary. Such a scenario would be permitted in model V, P and may be tested 
for by evaluating the number of posterior samples which satisfy b P > (1 — p); the definition of 
a grazing event. We find that (bp + p) < 1 to > 99.99% confidence and not a single posterior 
sample landed in this regime. Therefore, the hypothesis of an unblended grazing EB is highly 
improbable. 

Results: A Blended Grazing EB. A blended grazing eclipsing binary could feature in all 
four of the alternative hypotheses; V, C EB33M , "H^B33.6' ^bb 67.2 & ^eb 67.2- The first thing 
to note is that Bayesian evidence of all of these models is not significantly improved over the 
null hypothesis, Hp. In each of the four models, we checked the posterior samples satisfying 
a non-grazing configuration and found (bp + p) < 1 to > 99.99% in all four. Therefore, the 
hypothesis of a blended grazing EB is highly improbable. 

Results: A Blended 67.2 d EB. As discussed earlier, the 67.2 d period blended EB causes 
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a limb darkened transit but a flat-bottomed occultation. This difference in transit profile is ex- 
pected to make the scenario easier to distinguish than the 33.6 d blended EB. Indeed, this is what 
we found. Hypotheses V- c EBfi7 2 and "H|, B 67 2 yield A log Z = (— 22.7±0.5) and (— 133.1±0.5) 
respectively, relative to the null hypothesis Hp. This constitutes a QA-a and 16.1-er preference 
for the null hypothesis, for the two alternative hypotheses respectively. Therefore, the hypothe- 
sis of a blended 67.2 d EB is highly improbable. 

Results: A Blended 33.6 d Circular Orbit EB. With all previous scenarios now rejected, we 
are left with the blended 33.6 d EB only, which comes into two flavors: Heb^s & ^es,33.6- 
Let us here consider the circular case first. Strictly considering the Bayesian evidence results of 
our light curve profile analysis, there is no significant preference between the hypothesis of an 
unblended planetary transit and a blended 33.6 d EB. Note that these remaining blend scenarios 
include a Jupiter-sized planet transiting a blended main-sequence star. 

Scenarios involving a giant star or white dwarf can be excluded based upon the spectroscopic 
analysis discussed earlier. The ratio-of-radii in both V, C EB 33 6 is constrained to be p < 0.30 to 
> 99.99%. This means we must be dealing with a main-sequence star being transmitted by a 
smaller object, albeit with the possible presence of a blend. 

We discussed earlier how there is no detectable occultation in the data between the 33.6 d 
transits. Under the hypothesis of V, C EB ^ 3 6 , the other object must have (J^/J 7 *) < 0.0089 to 
3-er confidence. This indicates a small (p < 0.3), cool (T < 3000 K) object consistent with a 
planet or brown dwarf. 

Results: A Blended 33.6 d Eccentric Orbit EB. The final case of the eccentric 33.6 d EB 
yields a slightly improved Bayesian evidence over the planet-only model, with a significance 
of (2.4 ± 0.2)-cr. We do not consider this significant enough to overturn the planet hypoth- 
esis but nevertheless the scenario is considered plausible at this stage. A possible reason for 
the slight improvement is the ability of an eccentric orbit to generate a more diverse range of 
limb darkening profiles than the circular orbit case, when both models utilize quadratic limb 
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darkening. 

Results: Conclusions. We conclude that the only models which can adequately explain 
the spectroscopic analysis and the blend analysis are that of a planet transiting a star (Hp), a 
blended 33.6 d EB on a circular orbit (H% B ,33.e) an ^ a blended 33.6 d EB on an eccentric orbit 
(H e EB 33 6 ). In all cases a grazing transit configuration is excluded. 

A Dynamical Analysis The spectroscopic and blend analyses thus far leave us with three 
hypotheses. Here we will evaluate these hypotheses in light of the dynamical constraints from 
both the transit timing variations (TTV) and stability arguments. 

It has been established that the observed large transit timing variations cannot be caused by 
a moon. It is shown in the main text that the only plausible source for such large TTVs is a 
third body in the system orbiting close to the 5:3 orbital resonance. These fits do not allow us to 
directly measure the mass for the transiting body but do allow us to constrain the eccentricity to 
a high degree of confidence. Thus, the first hypothesis we may consider is the eccentric 33.6 d 
blended EB scenario, V, e EB 33 6 . 

Eccentricity Constraints. The spectroscopic analysis can be used to provide a mass and 
radius for the star using stellar evolution models, as was done earlier. However, if a substantial 
amount of blended light is present then the inferred properties would be unreliable. This is 
particularly salient in light of the faintness of the target and the subsequent lower-than-normal 
SNR spectra obtained. However, even in the case of a substantial amount of blended light, we 
are confident that the star is a dwarf on or near the main- sequence. We consider a wide range 
of corresponding plausible stellar masses to be 0.1 M < M * < 10.0 M Q . 

Using this mass range, we re-fit the TTVs each time varying all of the system parameters. 
After finding the maximum likelihood solution, we perturb the parameters in order to derive an 
upper limit on the orbital eccentricity of the transiting planet. Across the full stellar range, we 
find e P < 0.02 to 99.9% confidence. Note that the derived TTVs are insensitive to any blended 
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light and thus this limit is robust for B > 1. 

We conclude that the hypothesis of an eccentric 33.6 d blended EB is highly improbable. 
This now leaves only two remaining hypotheses to explain the data: Hp and 1~L%b 33 g. 

Mass Constraints. The two surviving hypotheses are identical in that they both have a body 
on a 33.6 d circular orbit eclipsing a main-sequence dwarf star on a non-grazing transit. The 
only difference is the amount of blended light (B = 1 versus B > 1). In both cases, the 
eclipsing object is small (p < 0.3 to 99.99% confidence) and cool (.Fp/J 7 * < 0.0089 to 3-a 
confidence, which is robust for any B value). The only plausible blend scenarios which could 
reproduce all of these constraints is a very cool M-dwarf, a brown dwarf or a larger planet. 
Clearly the planetary nature of the transiting body is not yet validated as these objects span a 
wide range of possible masses. 

The final and most powerful tool at our disposal is that the TTVs allow us to measure M c /M* 
for the third body. For a given M* then, two of the three masses in the system are known, along 
with their periods, semi-major axes, eccentricities and mutual inclination. The only unknown 
is the mass of the transiting object. We may place an upper limit on this value by iteratively 
increasing the mass until the system becomes dynamically unstable. 

Dynamical stability was investigated using the symplectic N-body code known as SyMBA 
(20), simulating the system for 1 Gyr using an integration timestep of 1.5 d. One can account for 
the possibility that the stellar mass derived from the spectroscopy is unreliable (due to blending) 
by investigating a wider, but plausible, range of stellar masses. To this end, we scanned the range 
0.8 M . < M. < 1.2 M . In order to be stable for 1 Gyr, we estimate \l t > < 5 Mj, <7Mj 
and < 9Mj for 0.8 M , 1.0 M and 1.2 M Q respectively, all of which exhibit a very sharp 
stability boundary. For our best-fit stellar mass of 0.9 M Q (see earlier spectroscopy discussion), 
the mass of KOI-872.01 is constrained to be M P < 6 Mj. 

We therefore conclude that the transiting object KOI-872.01 must be planetary in nature 
(since M P < 11 Mj\ (53)) and thus refer to the object as KOI-872b from here-on-in. KOI- 
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872c is validated as a planet based upon the precise measurement of its mass from the TTVs, 
specifically M c /M, = 3.97ig;^ x 1(T 4 , corresponding to M c = 0.376±g;^g Mj. 

A Search for Additional Transits Our analysis of the transit timing variations (TTVs) leads 
us to conclude that a second planet exists in the system with an orbital period of 57 d. Given 
that the inner planet transits, and conclusion of nearly coplanar orbits from the TTV fits, there 
seems a reasonable hope for detecting transits of a second planet. However, given the fact KOI- 
872c is Saturn-mass, it should be a gas giant with a sizeable transit, roughly of the same size as 
KOI-872b. Such a transit would be easily spotted even by eye but there is no evidence for such 
events. Fig. [S6] shows the phased data, accounting for the TTVs of KOI-872c upon the time of 
expected transit. Since the TTVs of KOI-872c depend upon the unknown mass of KOI-872b, 
we present 1 1 different realizations for various masses of planet b. In all cases no transit-like 
event is detectable. 

Despite the fact KOI-872c does not seem to transit, we initiated a search for additional 
transits in the system. The light curve was searched for transits using the Box Least-Squares 
method (54). After removing the transits of KOI-872b, we detected a significant signal (SNR~ 
14) in the light curve with an apparent depth of ~ 0.3 mmag, and a period of P = 6.7668293 d. 
The drop in brightness had a first-to-last-contact duration, relative to the total period, of q = 
0.0219, corresponding to a total duration of 3.6 hr. The transit candidate is hereby referred to as 
KOI-872.03. 

Transit Fit of KOI-872.03 Initial inspection of depth, duration and period for the KOI-872.03 
transits suggested a physically plausible signal. Given the high significance of the signal, we 
investigated further by performing a full transit light curve fit accounting for the limb darkening 
of the star, the variable diluted light factors and the finite integration time of the long-cadence 
data, all of which are ignored in the BLS search. 

We trim the data set to be within ±0.5 d of the expected transit times, as computed from 
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the BLS peak, in order to reduce the computation time. The 74 transit epochs require 74 OOT 
parameters to fit in conjunction with the transit parameters themselves. This large number of 
parameters makes a MultiNest fit unfeasible. Instead, we use a Markov Chain Monte Carlo 
(MCMC) routine with the Metropolis-Hastings rule. We executed two independent fits, one 
using a prior on p* from the Air fit and one with a uniform prior on p* between the boundaries 
for a main- sequence star. 

With the free p* fit, we detect a transit corresponding to a planet of size (Rp/R*) 2 = 
278^32 ppm (8.8 a). With the prior p* this becomes (R P /R*) 2 = 274+lppm (9.9 a), corre- 
sponding to a planet of size Rp = 1-701°;^ R 9 . The maximum likelihood realization of this 
latter fit is shown in Fig.[S7]and the corresponding parameters estimates are provided in Table 1 . 

As expected, the prior-p* fit retrieves virtually the same p* as derived from planet b, specif- 
ically p* = 1560^50 kg m~ 3 . The impact parameter converges to b = 0. 39^0^2 corresponding 
to i = (88.55t° il) . Curiously though, releasing this prior yields p* = 18201^ kg m~ 3 , 
b = 0.021q 44 and i = (88.92^11)° . In other words, without any prior information, the light 
curve yields a consistent stellar density (within 0.6 a), highly indicative that KOI-872.03 orbits 
the same star as KOI-872b with negligible eccentricity. 

Due to the low signal-to-noise, we were not able to determine individual transit times or 
durations for this object. We note that the expected TTVs of KOI-872b due to KOI-872.03 is less 
than 1 second in amplitude and thus undetectable with the current data. Further, the expected 
TTVs of KOI-872.03 are around 10 s, which are again too small to detect. Without TTVs, we 
cannot causally link the object to be transiting the same host star as KOI-872.03 (it could be 
transiting a background star). Even assuming it was in the same system stability limits allow 
the object to be as massive as 100 Mj and still be stable. Although formally unconfirmable, the 
derived p* suggests it is likely associated with the same star. 
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A Search for an Exomoon As discussed earlier, the observed TTVs cannot be adequately 
explained by an exomoon and only a second planet in the 5:3 resonance offers a valid solution. 
The presence of this second planet complicates our search for an exomoon. Further, the presence 
of stellar activity makes spot crossings and transit distortions probable, further exacerbating our 
search for a moon. Despite this, we here describe our efforts to search for an extrasolar moon. 

Fitting for a Moon Eclipse After Removing the TTVs, Mmti,mo- The first attempt we made 
was to forcibly remove the best-fit TTVs from the time series and then fit for a zero-mass moon. 
This fit allows for a finite radius moon and thus is merely a search for the moon eclipse. We 
perform two versions of the model fit; one setting the exomoon radius to zero and one allowing 
the parameters R s /Rp, (pp) 2//3 , Ps, 4>s, is an d to be freely varied (see (38) for various 
definitions of these terms). The two models are performed so that we have a null model to 
compare against. 

The results yield \ogZ = (12001.34 ± 0.37) for the null fit and \ogZ = (12025.89 ± 0.24) 
for the moon-transit fit, or a 6.7-er preference for the moon-transit model. Whilst certainly 
above our statistical significance threshold, one should recall that the star is active and these 
fitted events could merely be star spot crossings or activity -related events. Such events would 
be poorly sampled with the 30 minute cadence of the current observations though. 

The mass-ratio of the planet and the star may be determined using the light curve alone, as 
described in (55). For our estimate of the stellar mass, this allows us to compute M P directly. 
The results find that M P > 20 Mj for all modes, which exceeds the 6 Mj stability limit imposed 
on the system. We therefore conclude that none of these modes are genuine and most likely due 
to the presence of stellar activity. The likely presence of these spots also prevent us placing an 
upper excluded limit on a putative exomoon radius. 

Fitting for a Moon After Removing Planet TTVs, Mmt2,ro- We also tried removing the 
maximum likelihood TTVs from the planet-fit of KOI-872c. This is therefore a fit on the 
residual TTVs for a moon signal. Since the eclipse signal of the moon is likely unreliable 
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due to stellar activity, we limit this search to a model, where the moon is a point-mass (i.e. 
Rs = 0); model M.mt2,ro- We found this model unable to locate a significantly improved fit 
with A log Z = —(1.55 ± 0.44), relative to the null hypothesis. To illustrate this, Fig.lS8lshows 
the TTV residuals after removing the TTVs of KOI-872c along with the maximum likelihood 
model TTVs from M. mt2,ro- Although no exomoon is detected, we can use the results to place 
upper limits on a putative exomoon mass. This is particularly valuable given that radius limits 
are not possible due to the likely presence of spots. Fig. [S9] provides the corresponding mass 
limits, excluding Mg/Mp < 0.021 to 3-er confidence. 

TTV Constraints on Additional Planets When the computed TTVs corresponding to si 
are subtracted from the measured TTVs, this leaves a small residual signal with a ~1 minute 
amplitude. This is comparable to the measurement errors. 

The small amplitude of the TTV residuals can be used to place limits on the presence of 
additional planets in the system. Given that KOI-872b and KOT872c have nearly coplanar 
orbits, we tested a case in which the additional planet was placed in the invariant plane of the 
two confirmed planets. The orbits were followed by an iV-body integrator to see whether the 
computed TTVs are consistent with the residuals. 

The results are illustrated in Figs. IS 1 1 1 and IS 1 21 The small amplitude of residual TTVs 
provides an useful constraint on the third planet's mass and orbit. They rule out, for example, a 
Jupiter-mass planet on low-e orbit with 0.05 < a < 0.5 AU. 

Additional constraints can be obtained from the stability requirements. For example, a low- 
mass planet with e ~ and i ~ should have \a — a c \/a c > C(M C /M*) 2 / 7 , where M c and a c 
are the mass and semimajor axis of KOT872c, for the system to be stable (5(5), where C ~ 1.5 
(e.g., (57, 58)). For M c /M* = 4 x 10" 4 , this gives \a - a c \/a c > 0.160. 

A more accurate stability criterion was derived in (59). For M c /M* = 4 x 10~ 4 and e c = 
0.015, this criterion gives \a — a c \/a c > 1.8ey 5 (M c /M*) 1 / 5 = 0.162. The difference between 
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the two criteria is therefore negligible in our case. As the mean motion resonances become 
wider with planet's eccentricity, the eccentricity-dependent criterion of (59) should be used for 
larger e. 
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Figure SI: "Raw" (PA output) flux observed by Kepler from DR5 for Ql-6 of the source 
KOI-872 aka HCV/HCA-439. Overlaid is our model for the long-term trend, computed using a 
discrete cosine transform for each data set. Outliers and discontinuous systematic effects have 
been excluded. Transits (removed) marked with vertical gridlines. 
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Figure S2: Maximum likelihood model (black line) overlaid with the long-cadence Kepler data 
for KOI-872, surrounding the ramp affected transit. The simple exponential ramp model is 
shown to provide an excellent description of the instrumental effect. 
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Figure S3: x and y centroid positions of KOI-872 relative to the median value over a 0.8 day 
range, surrounding the transits of KOI-872b. Data temporally offset and phased to account for 
the transit timing variations. The vertical grid lines mark the first and fourth contact points. We 
find no deviation of the centroids between in- versus out-of-transit, which would have indicated 
a separated blend source. 
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Figure S4: By comparing the in-transit to out-of-transit centroid position for KOI-872b (solid) 
and KOI-872.03 (dashed), we find no evidence for a shift corresponding to a separated blend 
source. Our upper limits allow us exclude the blend factor, B, as a function of the separation of 
a companion. The top panel shows the results when plotted for B, the lower panel when plotted 
for magnitude difference. 
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Figure S5: YY-isochrone analysis of KOI-872. Using our spectrosopic observations and the 
constraint on p* from the transit light curve of KOI-872b, we are able to determine precise 
parameters for the host star. The backdrop of isochrones is for ages 0.2, 0.5, 1 .0, 2.0, . . . 13.0 Gyr 
(from bottom to top) and [Fe/H] = 0.41. The solution indicates an old star with 10 ± 3 Gyr age. 
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Figure S6: Using the TTVs of KOI-872b, we are able to predict the times of transit for KOI- 
872c. Assuming masses for planet b ranging from (top) to 6 Mj (bottom) for the primary in 
equal steps, we show the Kepler photometry phased upon these 1 1 candidate TTV ephemeres. 
Gray indicates the original data and black indicates 20-point binned data. We find no evidence 
for even a grazing geometry of KOI-872c. 
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Figure S7: Maximum likelihood realization of a Markov Chain Monte Carlo fit for the transits 
of KOI-872.03. Gray points show original data and black is 20-point phase binned data. The 
maximum likelihood model is in red and the residuals are shown below offset at +0.998. 
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Figure S8: Maximum a-posteriori moon fit to the TTVs of KOI-872b (using model Mmt2,ro), 
after removing the maximum likelihood TTVs due to the second planet, KOI-872c. We find no 
significant improvement by including a moon. 



51 



150 



CD 
"O 

_Q 

O 



100 




0.000 0.005 0.010 



0.015 0.020 0.025 0.030 

M S /M P 



Figure S9: Posterior distribution for the mass of an exomoon relative to the mass of KOI-872b 
as computed using model M. MT2,m, marginalized over the entire prior volume. We estimate a 
3-o- upper limit of M S /M P < 0.021 for this planet. 
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Figure S10: Confidence intervals. The dots denote the 99% confidence area around si (red 
triangles). We sampled the general neighborhood of si, determined x 2 for each parameter set 
(3 million parameter sets in total) and plotted a dot if x 2 < Xmm + A% 2 (99%) = 23.5. All 
parameters shown here are well constrained, including the orbital inclination of KOI-872c (the 
gray area in (d) is ruled out by the lack of KOI-872c's transits). This result was obtained while 
fixing e h = 0. A nearly identical result was obtained by letting e h (and w b ) vary. In that case, 
X 2 < Xmin + Ax 2 (99%) with Ax 2 (99%) = 16.9 for 15 - 9 = 6 DOF gives e b < 0.02, leaving 
Wb unconstrained. The scaled semimajor axis in panel (b) is defined as (a — a c )/a c , where a c is 
the best-fit semimajor axis value for solution si. 
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Semimajor Axis (AU) 

Figure SI 1: Constraints on the presence of an additional Jupiter-mass planet. The dashed area 
shows where the Jupiter-mass planet would induce the TTV amplitude of KOI-872b in excess 
of 1 minute. If the residual TTVs are caused by a Jupiter-mass planet, this planet should have a 
and e near the boundary of the shaded area. 
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Figure S 12: Constraints on the presence of an additional planet on circular orbit. The line shows 
where the planet induces the TTV amplitude of 1 minute. Note that the low-mass planet on a 
circular orbit exterior to KOI-872b cannot explain the TTV residuals, because such a planet 
would be perturbed by KOI-872c and would be dynamically unstable. 
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Table SI: Priors used for various free parameters in the MultiNest fits. U[x, y] is a uniform 
prior between x and y. Af[x, y] is a Gaussian prior with mean x and standard deviation y. S[x] 
is a delta-function prior centered on x. J'[x, y] is a modified Jeffrey's prior with an inflection 
point at x and a maximum limit of y. We use the replacements Tq = 2455053.2815 BJDutc. 
Pp = 33.6013 d and r* = r * + nPp. a photo is the median photometric error outputted from the 
Kepler pipeline for the time series under analysis. * = for Mv fits these terms are independently 
fitted to each transit epoch. N/A = parameter is fixed to some arbitrary value, since it has no 
influence on the fits, t = for Mm,ro an d M.mt2,ro, this term is replaced with (M S /M P ) and 
treated as uniform between and 1. x = for M.mti,mo, this term is fixed to zero to represent a 
zero-mass moon. 
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All TTVs removed 
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Planet TTVs removed 


•M.MT2,iml\ 


(12006.43 ±0.37) 
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Table S2: Bayesian evidences for the planetary system family of models fitted to the KOI- 
872 photometry. The data favor model Mr, a planet-only model with transit timing variations 
(TTV). 
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Table S3: Epoch-specific fitted parameters for the KOI-872 system. Results computed from 
the weighted posteriors resulting from model Mr, using MultiNest. TTVs relative to 
maximum a-posteriori linear ephemeris derived in model fit Mp\ Pp = 33.6013506 d and 
r = 2455053.2815010. Physical system parameters are provided in Table 1. 
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Table S4: Epoch-specific fitted parameters for the KOI-872 system. Results computed from the 
weighted posteriors resulting from model My, using MultiNest. TDVs relative to median 
duration of all duration realizations in M v ; T 14 = 250.5 mins. Physical system parameters are 
provided in Table 1 . 



Hypothesis 



Z 



a Preference over "Hp Blend analysis conclusion 



■Hp 

njc 

n EB,33.6 
nje 

rl EB,33.6 
njc 

n EB,67.2 
nje 

n EB,67.2 



(12000.13 ±0.37) - Plausible 

(12000.95 ±0.37) (0.8 ± 0.4) Plausible 

(12004.29 ±0.35) ±(2.4 ±0.2) Plausible 

(11977.45 ±0.39) -(6.41 ±0.08) Excluded 

(11867.05 ±0.36) -(16.13 ± 0.03) Excluded 



Table S5: Bayesian evidences for a blend family of hypotheses fitted to the KOI-872 photom- 
etry. Out of the various blend scenarios tried, we find the light curve shape has no significant 
preference between hypotheses Hp, H C EB336 and H e EB 33 6 . 
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